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Abstract 

Parton distribution functions (PDFs), which describe probabihty densities of quarks and gluons in 
the proton, are essential for interpreting the data from high energy colliders. The PDFs are measured by 
approximating them with functional forms that contain many adjustable parameters. Those parameters 
are determined by fitting a wide variety of experimental data. This paper examines the uncertainty 
that arises in this procedure from the choice of parametrization, and shows how that uncertainty can 
be reduced by a technique based on Chebyshev polynomials. 
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1. INTRODUCTION 



Interpreting the data from high energy colhders such as the Tevatron and LHC requires parton 
distribution functions (PDFs), which describe the probabihty density for quarks and gluons in 
the proton, as a function of hght-cone momentum fraction x and QCD factorization scale /i. In 
current practice jl, 2|, the PDFs are measured by parametrizing them at a low scale /xq, using 
functional forms in x that contain many adjustable parameters. The distributions at higher fi 
are calculated using QCD renormalization group equations, and the best-fit parameter values 
are found by a "global analysis," in which data from a variety of experiments are simultaneously 
fitted by minimizing a measure of fit quality. 

In the Hessian method [sl, the uncertainty range of the PDFs is estimated by accepting all fits 
for which is not more than some fixed constant Ax^ above the best-fit value. Traditionally, 
AX^ . 100 a or Ax^ ^ 50 have been used. When these "large- A^^ values are used, weight 
factors or penalties are included to maintain an adequate fit to each individual data set py]; or 
the uncertainty range along each eigenvector direction can instead be estimated as the region 



where the fit to every experiment is acceptable 



?he large Ax tolerance criterion has long 



6|, Q] eschew it in favor of the Ax^ = 1 for 



been a source of controversy, and some PDF groups 
68% confidence (Ax^ = 2.7 for 90% confidence) that would be expected from Gaussian statistics. 

A potential motivation for large Ax^ is based on conflicts among the input data sets, which 
signal unknown systematic errors in the experiments, or important theoretical errors introduced 
by reliance on leading-twist NLO perturbation theory. It would make sense to scale up the 
experimental errors to allow for such conflicts Q]. However, the conflicts among experiments 
were recently shown QQ to be taM. sn,a,l: .ea..ed d.c„ betweeu the ..phcat.o. 
from each experiment with respect to all of the others suggest a minimum A^^ ~ 10 for 90% 
confidence, but supply no incentive for Ax^ 100. 

Another source of uncertainty in PDF determination is the parametrization dependence error, 
which comes from representing PDFs at /iq — unknown continuous functions — by expressions that 
are adjustable only through a finite number of free parameters. In traditional practice, one adds 
flexibility to the parametrizing functions one parameter at a time, until the resulting minimum 
X^ ceases to decrease "significantly." However, wherever one chooses to stop adding fitting 
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parameters, a further small decreases in remain possible. This aspect of the PDF problem, 
namely that the number of fitting parameters is not uniquely defined, spoils the normal rules 
such as Ax^ = 1 for 68% confidence, which would otherwise follow from standard Gaussian 
statistics. This point is illustrated in Sect. [2] by some hypothetical examples. A method that 
invokes Chebyshev polynomials to increase the freedom of the parametrization is introduced in 
Sect. 131 and applied to a typical PDF fit in Sect. IH Conclusions are presented in Sect. [51 



2. HYPOTHETICAL EXAMPLES 



Let z be the displacement from the minimum point in along some specific direction in 
parameter space. It can be normalized such that 



in the neighborhood of the minimum. 



X=z'^C (1) 
The parameter z can be any one of the eigenvector 

CI □ 

coefficients Zi that are discussed in [9( or [10|. Or by means of a suitable linear transformation 
X = a + f5z^ z can represent the prediction for some cross section X that depends on the PDFs; 
or simply a PDF itself for some specific flavor, x, and ^. According to standard statistics, 
Eq. dll) would imply z = ± 1 at 68% confidence and ± 1.64 at 90% confidence. If we assume 



instead — guided by [lO| — that Ax^ = 10 for 90% confidence, we would expect z = ± 3.16 
at that confidence. However, the following argument shows that the uncertainty range could 
actually be much broader than that. 

Suppose that, in order to reduce the dependence on the choice of parametrization, we in- 
troduce additional flexibility into the PDF model through a new parameter y, which is defined 
such that x^ reduces to Eq. at y = 0. (To achieve a substantially improved fit, it will likely 
be necessary to increase the flexibility in more than one flavor, and therefore it will be necessary 
to introduce several new fitting parameters. The parameter y therefore represents displacement 
in a direction that is defined by some particular linear combination of new and old parameters.) 

Figure [H shows two hypothetical contour plots for x^ as a function of y and z. The contour 
interval is 10. In each case, introducing the new parameter reveals that 2 = 10 is a better 
estimate of the true value of z, so the prediction according to y = 0, that 2; = ± 3.16 at 90% 
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FIG. 1: Contours of = 3010, 3020, 3030, ... in two hypothetical examples. In each case, the best fit 
has = 3000 at z = 10. Meanwhile, if the fit is confined to y = 0, the Ax^ = 10 error limits appear 
to be z = lb 3.16, which is far too restrictive. 

confidence, is incorrect. In the scenario of the left panel, the additional freedom measured by y 
has reduced the best-fit x"^ by 50; while in the right panel, the reduction is only 5 — a change so 
small that one might easily have been content to mistakenly settle for y = 0. 

In the hypothetical examples of Fig. [H Ax^ = 10 yields an estimate of uncertainty for fits 
with y = that is far too narrow — even though in one case, the additional freedom only allows 

to be lowered by 5. Appendix I shows that the qualitative form of the dependence of on 
y and z shown in Fig. [1] arises rather generally, when additional freedom is introduced into the 
par ametrizat ions. However, it remains to be seen whether such large quantitative changes arise 
in typical PDF fitting. A new parametrization method introduced in the next Section will be 
used to answer that question. 
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3. CHEBYSHEV PARAMETRIZATIONS 

In a recent typical PDF fit (CT09 jit), the gluon distribution at jjo was parametrized by 

= aox"i (l-x)'^^eP(^) (2) 

where 

p{x) = Osa/x + a^x + a^x^ . (3) 

The same form was used — with different parameters of course — for Uy = u — u and = d — d, 
except that 03 was set to in because that distribution is less constrained by data. 

To provide greater flexibility in the parametrization, it would be natural to replace p{x) by 
a general polynomial in i/a: : 

n 

p{x) = E b, x^/' . (4) 
i=i 

This form has several attractive features: 

1. The power-law dependence in x at x — )■ 0, with subleading terms suppressed by additional 
powers of approximately x°'^, is expected from Regge theory. 

2. The power-law suppression in (1 — x) at x — )■ 1 is expected from spectator counting 
arguments. 

3. The exponential form e^'-^-' allows for the possibility of a large ratio between the coefficients 
of the power-law behaviors at x — )■ and x — 1, without requiring large coefficients. It 
also conveniently guarantees that g{x) is positive definite — although that could in principle 
be an unnecessarily strong assumption, since the MS parton distributions are not directly 
observable, so it is only required that predictions for all possible cross sections be positive. 

4. Restricting the order n of the polynomial in Eq. (jl]) can express the assumed smoothness 
of the parton distributions — although if n is large, some additional condition must still be 
imposed to prevent rapid variation. 

The constraints on smoothness and behavior at x — and x — 1 are important. Without 
them, the momentum sum rule and the valence quark number sum rules would have no power, 
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because mildly singular contributions near x = and x = 1 could make arbitrary contributions 
to those integrals over x, without otherwise affecting any predictions. 

In past practice, only a small number of non-zero parameters bj have been retained in 
as exemplified by the typical choice (|3]). The number of parameters can be increased to add 
flexibility, and thereby reduce the dependence on the choice of form for parametrization. How- 
ever, that quickly runs into a technical difficulty: as more fitting parameters are included, the 
numerical procedure to find the minimum of becomes unstable, with large coefficients and 
strong cancellations arising in p{x). The resulting best fits, if they can be found at all, contain 
implausibly rapid variations. 

This technical difficulty can be overcome using a method based on Chebyshev polynomials. 
These polynomials have a l ong tradition in numerical analysis, but have only recently begun to 



be apphed in PDF studies They are defined — and conveniently calculated — by recursion: 

To{y) = 1, T,{y) = y 
Tn+i{y) = '^yTniv) - Tn-i{y) . (5) 

Since Tj(y) is a polynomial of order j in y, the parametrization (jl]) can be rewritten as 

n 

p{x) = Y,CjTj{y) , (6) 

where y = 1 — l^fx conveniently maps the physical region 0<x<lto —\<y<\. 

The parameters Ci, . . . , c„ are formally equivalent to the parameters 61, . . . , but they are 
more convenient for fitting, because the requirement for smoothness in the input PDFs can be 
expressed by forcing the Cj parameters to be reasonably small, particularly at large order. This 
can be seen from the following property of the Chebyshev polynomials: 

Tj{y) = cos{j9) where y = cos9. (7) 

With the mapping y = 1 — l^fx (or more generally y = 1 — 2x°'), Tj{y) has extreme values of ±1 
at the endpoints and at j — 1 points in the interior of the physical region < x < 1. Chebyshev 
polynomials of increasingly large j thus model structure at an increasingly fine scale in x. A 
particular way to force the coefficients cj to be small, especially at large j, in order to enforce 
smoothness, is described in Appendix II. 
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Using the Chebyshev polynomial method, it becomes possible to obtain fits, and even eigen- 
vector uncertainty sets, with two to three times more free parameters than were tractable in 
previous PDF fitting. The method is applied in the next Section to evaluate the parametrization 
error in a traditional fit. 



4. FITS USING THE CHEBYSHEV METHOD 




FIG. 2: Uncertainty of gluon distribution, and d and u quark distributions, calculated at Ax^ = 10 
using the CTEQ6.6 form of parametrization with 22 parameters, compared to fits using the Chebyshev 
polynomial method with 26 parameters (long dashed), 31 parameters (medium dashed), and 36 pa- 
rameters (short dashed). The two versions of each curve (with or without dots) correspond to different 
choices for the smoothness penalty. 



Figure [2] shows the fractional PDF uncertainty obtained using the parametrization method of 
CTEQ6.6 which has 22 free parameters. The uncertainty limit is defined here by Ax^ = 10, 
which is the range suggested by observed conflicts among the input data sets lOj. The curves 
show the ratio to that central fit, of fits made using the Chebyshev polynomial method defined 
specifically in Appendix II, with 26, 31, and 36 free parameters. Two fits are shown for each of 
these choices, corresponding to different strengths of the penalty that is added to to enforce 
smoothness. Introducing more flexible parametrizations is seen to shift the best-flt estimate of 
the PDFs by an amount that is in some places comparable to or somewhat outside the range of 



7 



uncertainty that was estimated on the basis of conflicts among the input data sets. 




FIG. 3: Same as Fig. [2l but using a linear scale in x to show the behavior at large x. 



Figure |3] displays the same results over a wider range in the ratio, using a linear scale to reveal 
the behavior at large x. In the region x > 0.8, where the PDFs are very poorly constrained, the 
improved parametrizations sometimes lead to flts that lie far outside the original uncertainty 
band. This comes about because the uncertainty range was obtained using the Hessian eigen- 
vector method, which is based on assuming quadratic behavior (second-order Taylor series) of 

as a function of the fltting parameters. For quantities that are weakly constrained (i.e., along 
"flat directions" in parameter space), this approximation can easily break down before has 
increased very much. 

Table [T] lists fo^^ the flts shown in Figs. [2] and [31 All of these flts would be acceptable 
according to a crude "hypothesis testing" criterion, since is close to the total number of 
data points (minus the number of fltting parameters) — well within the range ± y^2iVp^ . The 
additional parametrization freedom allowed by the Chebyshev method produces a distinctly 
lower overall x^y ^ind it therefore can be assumed to provide a global best flt that is a better 
approximation to the true PDFs. This is analogous to moving from the minimum along 
y = to the minimum in the {z, y) plane depicted in Fig. [1] Using more than 36 parameters 
for this data set offers negligible further improvement in the flt, since beyond that number, the 
further decrease in x^ becomes smaller than the number of added parameters, so x^ per degree 
of freedom no longer decreases. 
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Parametrization 


Npar 


Penalty 




CTEQ6.6 form 


22 





3085.9 


Chebyshev 


26 


2.6 


3083.0 


Chebyshev 


26 


0.2 


3082.7 


Chebyshev 


31 


5.2 


3066.5 


Chebyshev 


31 


1.2 


3061.8 


Chebyshev 


36 


6.6 


3055.3 


Chebyshev 


36 


0.8 


3051.3 



TABLE I: Number of fitting parameters Npar, smoothness penalty added to x^? ™d for the fits 
shown in Figs. [2] and [3] to 3027 data points. 

5. CONCLUSION 



Improving the flexibility of the input parametrizations used in PDF analysis, using a method 
based on Chebyshev polynomials, reveals a significant uncertainty in previous PDF determi- 
nations caused by parametrization dependence. In regions where the parton distributions are 
fairly accurately constrained by experiment, the uncertainty due to parametrization is roughly 
comparable to the uncertainty suggested by observed conflicts between different data sets that 
are used in the global analysis. In regions where the experimental constraints are weak, such as 
the large-a: regions shown in Fig. [3l the parametrization-dependence uncertainty is considerably 
larger. This comes about because the uncertainty there comes mainly from parametrization 
dependence, and the reliance of the Hessian method on the quadratic approximation naturally 
underestimates such uncertainties. 

The parametrization uncertainty must be combined with the uncertainty evidenced by con- 
flicts among the input data sets. The latter has been previously estimated to be given by 



Ax^ = 10 at 90% confidence 10|. The parametrization uncertainty found here is on the same or- 
der, so combining the two uncertainties naively leads to a prescription of Ax^ = 20 or Ax^ = 40, 
depending on whether the uncertainties are added in quadrature or linearly. Because there are 
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additional uncertainties that have not been treated here, such as the choice of the effective mass 
of the charm quark, and the choice of which data sets to include in the analysis, the use of Ax^ 
in the range 50 — 100 for 90% confidence continues to appear reasonable. This estimate is also 
supported by recent fitting, which will be reported elsewhere, in which the many DIS data sets 
from the HI and ZEUS experiments at HERA are replaced by a single joint data set [3] in 
which the experimentalists have reduced the systematic errors by combining results from the 
two experiments. The best fit using the new data set is different from the best fit using the old 
data by an amount corresponding to A^^ = 55. 

Further confirmation of the size of the uncertainties from a practical point of view is shown 
in Fig. m which compares fits shown in Fig. [2] with the original CTEQ6.6 j4| for d{x) at /i = 
100 GeV. The wide band is the CTEQ6.6 uncertainty (weighted Ax^ = 100). The narrow band 
and solid curve are the new fit using the CTEQ6.6 parametrization (22 parameters). The central 
fit differs from CTEQ6.6 because of additions to the input data (inclusive jet and Z° rapidity 
distribution data from run II at CDF and DO), the introduction of a mild cut, {l—x)Q^ > 9 GeV^, 
on fixed-target muon DIS data as suggested by results of [l^, and because all experiments are 
assigned weight 1 in the new fit. The width of the uncertainty is smaller, mainly because 
Ax^ = 10 was assumed here. The dot-dashed curve is the most extreme of the Chebyshev fits 
shown in Fig. [2] (36 parameters). Note that the difference between this curve and the result 
of the CTEQ6.6 parametrization is nearly as large as the uncertainty originally estimated in 
CTEQ6.6. This proves that the CTEQ6.6 parametrization form induces an uncertainty due to 
parametrization at the level of Ax^ = 100. It is reasonable to expect that the parametrization 
uncertainty is somewhat smaller than that in the case of the Chebyshev fits, since they are more 
flexible. 



An alternative to the Hessian method is provided by the NNPDF approach 13|, which 
represents the PDFs at fiQ by a Neural Network model that contains a very large number of 
effective parameters, and thereby avoids parametrization dependence almost entirely. Broadly 
speaking, the uncertainties estimated by NNPDF are similar in magnitude to the uncertainties 
estimated previously by the Hessian method using A^^ ~ 100. A more detailed comparison, 
which requires a discussion of differences in assumptions based on nonperturbative hadronic 
physics — such as positivity and behaviors in x — ?■ and x — )■ 1 limits — will be undertaken in a 
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X 



FIG. 4: Uncertainty of d at /x = 100 GeV, relative to CTEQ6.6 [J]. Wide band is the CTEQ6.6 uncer- 
tainty (weighted A^^ = 100). Narrow band and sohd curve are the fit using CTEQ6.6 parametrization 
(22 parameters, Ax^ = 10) shown in Fig. [2j Dot-dashed curve is a fit with 36 parameters, also shown 
in Fig. [21 



future work. 

The sensitivity to parametrization choice demonstrated in Figs. [2] and [3] shows once again 
that uncertainty estimates based on the classical Gaussian "Ax^ = 1" estimate, with a fixed 
choice of parametrization, are overly optimistic — especially since the groups who apply that 
criterion generally use PDF parametrizations with substantially less fiexibility than any of those 
used here. (In some recent studies, the parametrization dependence has been incorporated in 
a simple way, by expanding the estimated uncertainty r ang e to include results from each of a 
variety of specific choices for the parametrization forms 14j.) 

The parametrization effects found here are comparable to, or slightly larger than, the uncer- 
tainty estimated by a Ax^ = 10 analysis. They are not as large as the hypothetical examples 
whose spectre was raised in Section O except perhaps in regions where the uncertainties are al- 
ready very large. In future analyses, the Chebyshev method may be combined with the Hessian 
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method, to compute PDF uncertainties that are less subject to parametrization dependence. 
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Appendix I: General form of near its minimum 

This Appendix shows that the qualitative behavior of y^iy^z) hypothesized in Fig. [1] arises 
under rather general assumptions. 

Let us assume as usual that can be approximated in the neighborhood of its minimum 
by Taylor series through second order. To find the uncertainty of a particular variable, we can 
assume by means of a linear transformation of that variable, that we are interested in the value 
of a parameter z for which = + C at ?/ = 0, as in Eq. ([T]). Now let y represent an additional 
fitting parameter, which was previously held fixed at 0. By Taylor series, the expression for 
expands to become 

= ^2 ^ ^2 ^ 2A^y + 2By + C, (8) 

where the coefficient of was chosen to be 1 without loss of generality, by scaling that variable. 
Eq. (IH]) implies that the contours of constant are ellipses whose major and minor axes make 
an angle of ±45° with respect to the y and z axes, as in the specific examples of Fig. [TJ The 
ratio of minor axis to major axis of the ellipse is — |A|)/(1 + \A\) ; and \A\ < 1 is required, 
since x^ must have a minimum. The minimum of x^ occurs at 

^0 = - A^) , = -5/(1 - , (9) 

and its value there is 

= C - D, where D = 3^(1 - A^) . (10) 
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Relative to the y = situation given by Eq. ([T]), introducing the additional parameter y thus 
allows the best-fit to be lowered hj D. At the same time, it demands that the uncertainty 
range for z be extended at least far enough to include Zq, and hence demands Ax^ > Zq in the 
y = model. 

The hypothetical examples shown in Fig. [T] correspond to A = 5/5 = — ^2/3 =^ {zq = 10, 
D = 50), and A = 2B = - ^20/21 =^ {zq = 10, D = 5). Whether or not such large shifts 
occur in typical PDF fitting is the subject of Sect. HI 

Appendix II: Details of the Chebyshev method 

This Appendix describes a specific method that uses Chebyshev polynomials to parametrize 
the PDFs in a manner that allows considerable freedom in the parametrization, while requiring 
the functions to be smooth in a reasonable sense. 

Each flavor is parametrized using the form ([2]) at scale /io, with p{x) given by ([6]) using 
y = 1 — I^Jx. The Chebyshev polynomial coefficients Cj are scaled in a manner such as 

c^ = d,l3\ (11) 

to define fitting parameters rfj for which ^ represents an appropriate condition of 

smoothness. This smoothness condition is enforced in the fitting procedure by adding a penalty 
k ^d'j to x^i with the coefficient k chosen to make the total penalty on the order of 1. These 
penalties are listed in Table [H 
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